生存分析操作指南, 从医学中出来的社科瑰宝

可有偿投稿计量经济圈，计量相关则可

邮箱：[e-mail]

所有计量经济圈方法论丛的do文件, 微观数据库和各种软件都放在社群里.欢迎到计量经济圈社群交流访问.作者:@alphabeta

推荐阅读：

3.计量经济圈经济社会等数据库合集, 社科研究的米

今天，我们“面板数据研究小组”将为圈友们引荐关于“生存分析”(survival analysis)的操作和运用。若需要完整do file和数据，请到计量经济圈社群提取使用。

生存分析是研究生存现象和响应时间数据及其统计规律的一门学科。通常生存数据都具有一个重要的特点: 在研究期间结束时，某些个体身上还未出现令人关心的事件，即存在删失(censored)的现象。比如我们想要研究1990-2003年期间某些企业的存续时间，可在此时间范围之外这些企业的生存状况我们并不知道。如果这些企业成立在1990年以前，我们对企业成立到1990年这一时间段的生存状况就无从考察，而对这一问题进行忽略，就会导致数据删失问题，从而造成对样本企业生存时间错估。在这种情况下，研究者通常是将成立在1998年之前的企业数据删除，只保留成立在1990-2003年期间的企业，这样就去掉了左删失的观测样本。而对于 2003年仍然存在的这些企业，我们同样不能预测其之后的生存状况，即数据的右删失问题，这个可没法删除(数据都没有咋删除)。

对该类问题的分析并不能简单套用通常的统计方法，不过生存分析通常可以很好地解决这一“右删失”问题。在Cox比例风险模型中，我们定义企业的生存为企业成立到退出市场的状态(即是个二值变量)，企业的生存时间(survival time)定义为企业成立到退出市场所经历的时间(即是个连续或离散变量)， 将企业退出市场的事件称之为“失败”(failure)。一个企业生存时间的分布可以通过以下两种方法来描述: 一是反映个体存活超过时间ｔ的概率的生存函数(survival function)，后面会通过K-M曲线展示；二是反映在t-1时刻的个体下一瞬间经历某特定事件概率的危险率函数(hazard function)，实际上也有类似K-M的退出曲线，不过下文没有做进一步的展示。

在多因素生存分析中目前应用最为成熟的是Cox比例风险模型 (proportioanl hazard model)—Cox PH模型，它能有效地解决样本数据的右删失问题，准确地反映各影响因素与生存时间之间的关系。下面，我们就首先使用半参数估计Cox-PH模型开展生存分析，后面再用参数模型AFT做一个类似的生存分析，看看两者在多大程度上可以互相支持一个结论。

首先，在生存分析之前需要做的是定义生存分析数据(使用的程序是stset)。稍微讲解一下stset下面的各个选项。下面的例子使用的数据是多期面板数据，因此需要使用id来标明使用的各个观测值的id。time就是咱们的观测(记录)时间或年份。failure(m=1)表明，咱们这些个体在m=1的时候进入“failure”状态(在m=others值的时候就是一种删失状态)。origin()表明，这些个体在什么时候开始进入危险状态。这里的意思就是，从进入危险状态到最后死亡或退出，毕竟还有一个需要经历的过程。

stset time, id(id) failure(m==1) origin(stime) //进行了生存数据的定义

下面，我们绘制一个K-M生存曲线，该图表示观测对象生存时间大于时间t的概率。从下图可知，随着分析时间的变长，该个体生存概率实际上是在降低，退出或者死亡的可能性大大提高了。个体的生存函数在其生存时间为1至4年之间曲线下降最快，即该个体的生存时间主要集中在1到4年。

我们还可以绘制二条生存曲线，来大致看一看比如个体A与个体B的生存分布是否存在显著差异。在经济学中，我们可以这样设想，东部企业与西部企业的生存分布是否存在显著差异。下方的生存曲线图显示，这两者确实存在差异，当e=1的时候的生存分布曲线始终要比e=1的时候的生存分布曲线低。但，是否两者存在统计上的显著性，需要通过如下的一些检验方法。

比较生存曲线常用的检验方法有对数秩检验(Log-rank test)、威尔克森检验(Wilconxon test)和似然比检验(Likelihood ratio test)。当生存时间的分布为威布尔分布或属于比例风险模型时，对数秩检验(Log-rank test)效率较高；当生存时间的分布为近似对数正态分布时，威尔克森检验(Wilconxon test)效率较高；当生存时间的分布近似呈指数分布时，似然比检验(Likelihood ratio test)效率较高。

下面的Log-rank test检验显示，e=1与e=0两者的生存分布存在显著的差异。在实证研究中，你可以检验那些有阿里巴巴背景的被投资企业与有腾讯背景的被投资企业是否在生存概率方面存在显著差异。

下面是使用Cox-PH模型做回归所得到的结果。因为一些原因，这里的y是核心解释变量，千万不要把他当成是因变量(给你造成的不便，表示歉意)。下面显示的系数是hazard ratio，我们最好通过换算得到最终的系数并进行解释。

这是咱们通过hazard ratio换算得到的系数，从中可以看到，核心解释变量y的符号为正，因此通过提高y的水平可以让这个个体的风险水平升高，即生存的可能性更低。

但是运用Cox比例风险模型需要数据样本满足两个基本假定：一是等比例风险假定，任意个体风险函数之比随时间改变；二是对数线性假定，即对于模型中的连续变量，任一个体的对数风险与协变量呈线性关系。

在实证研究中，第一个等比例风险的假定是必须做的，这个与咱们之前学习的ordinal logit有点相似之处，因为一旦违反假定就会导致估计结果出现偏差。下面做的是等比例风险假定，值得我们关注还是最后一列各个协变量的P值和整体P值。除了x这个协变量外，包括咱们核心解释变量y在内的所有其他协变量都是满足等比例风险假定的。而且，整体的P值也显示这个模型是满足等比例风险假定的，因此前面的结论是可信的。

接着，我们加入一个新的选项——enter()，这个选项与前面的failure()选项是对应的，表示我们从m=1开始到m=2这个阶段死亡或衰退。origin在这里与上面的意思一样，代表在处于危险的时间。

现在，我们看看下面的回归结果(注: 这里直接得出的是系数，而不是hazard ratio)。核心解释变量y不显著了，表明在从m=1到m=2这个阶段，y并不能让这个个体的风险水平降低，也不能延缓其衰退速度。

下面的结果表明该生存模型满足等比例风险假定，因此该模型本身是没什么问题，看来确实是因为y在m=1到m=2这个区间段影响不了该个体的生存概率。

在多因素生存分析中，COX-PH比例风险模型和加速失效时间模型(AFT model)是较为常用的两个模型，均能有效地解决样本数据的右删失问题，从而准确地估计出生存时间与各影响因素之间的关系。但是使用 COX 比例风险模型需要数据满足一定的假设，要确保任意个体风险函数之比不随时间改变，即前面的等比例风险假定。在前面的表格中，我们通过等比例风险检验发现，x这个协变量并未通过该“假定”，因此是违反了COX-PH模型的assumption。在这种情况下，我们也可以考虑使用一种参数化的生存分析模型——加速失效时间模型。

在加速失效时间模型中，存活时间的自然对数logT，可由一组协变量的线性函数表达出来，进而构造出一个线性模型: logTj = Xj*β + Zj 。其中Xj 是一组协变量，β是这组协变量的回归系数，Zj表示密度函数为f(·)的误差。误差项不同的分布决定了不同的回归模型，如果f(·)是正态分布( normal)密度函数，那么该模型就是lognormal的回归模型; 同理，如果f(·)是logistic 密度函数，该模型则为loglogistic 回归模型; 如果f(·)是极值密度函数，该模型则为exponential回归模型或者Weibull回归模型 。

加速失效模型的效果取决于 exp(-Xj*β) 对于时间范围改变的程度，当exp(-Xj*β) 大于1时，时间就是加速的; 当exp( -Xj*β) 小于 1 时，时间就是减速的。也就是说，当观测个体在 t 时刻基准的生存概率为S( t) ，那么具有 Xj 协变量的观测个体在 t 时刻的生存概率S(·) 应该等于在exp(-Xj*β) t 时刻的估计值。所以，当 Xj*β 增加时，并不意味着时间的加速，相反是时间的减速，或者增加了等待失效( 失败) 的期望时间，也就是生存时间会增加。

在AFT模型中，我们可以通过log似然值和AIC, BIC等信息标准去比较哪种误差分布形态更加好地拟合了该数据。与平时我们使用的标准一样，log似然值越大越好，而AIC, BIC值越小越好。

下面就开始做AFT模型的回归。第一步还是需要定义咱们的数据是生存数据，与前面的一样。

stset time, id(id) failure(n==1) origin(stime) 

log normal AFT回归的结果如下，核心解释变量在这里并不显著。

Weibull AFT回归的结果如下，核心解释变量在这里也并不显著。

从上面的Cox-PH模型与AFT模型，我们发现两者还是有很大的区别。在Cox-PH中显著的核心解释变量y，却在AFT模型中变得不在显著，这至少证明咱们的稳健性还是有待商榷。

时间有限，就先引荐到这里，若对完整do file和数据有需求，可以到社群提取。

所有计量经济圈方法论丛的do文件, 微观数据库和各种软件都放在社群里.欢迎到计量经济圈社群交流访问.

可以到计量经济圈社群进一步访问交流各种学术问题，这年头，我们不能强调一个人的英雄主义，需要多多汲取他人的经验教训来让自己少走弯路。

计量经济圈当前有几个阵地，他们分别是如下4个matrix：

①小鹅社群：数据软件书籍等所有资料(最多)，

②微信群：服务于计量经济圈社群群友(最活跃)，

③研究小组：因果推断, 空间计量, 面板数据(最专业)，

④QQ群：2000人大群服务于社群群友(最大)。

计量经济圈是中国计量第一大社区，我们致力于推动中国计量理论和实证技能的提升，圈子以海内外高校研究生和教师为主。计量经济圈六多精神：计量资料多，社会科学数据多，科研牛人多，名校人物多，热情互助多，前沿趋势多。如果你热爱计量并希望长见识，那欢迎你加入到咱们这个大家庭(戳这里)，要不然你只能去其他那些Open access圈子了。注意：进去之后一定要看小鹅社群“群公告”，不然接收不了群息，也不知道怎么进入咱们独一无二的微信群和QQ群。

只有进去之后才能够看见这个群公告